The Need for Spatial Analysis

Spatial statistics asks ‘How much is where’ instead of the traditional ‘How much?’ The key element to the analysis is the autocorrelation of the observations in space.

Autocorrelation: correlation between elements of data separated by a given interval of space/time (auto = correlation of a variable with itself)

Correlation: interdependence of variable quantities (traditional statistics, longitudinal data, repeated measures)

Types of Spatial Data

(Classification discussed here according to Cressie (1993))

Consider a spatial process in \(d\) dimensions: \(\{Z(s): e \in D \subset \mathbb{R}^d\}\). \(Z\) is the attribute we observe at spatial location \(s\) (a vector if \(d\) co-ordinates). Most often \(d=2\) - Cartesian co-ordiantes.

library(spatstat)

## Loading required package: spatstat.data

## Loading required package: nlme

## Loading required package: rpart

## 
## spatstat 1.59-0       (nickname: 'J'ai omis les oeufs de caille') 
## For an introduction to spatstat, type 'beginner'

a <- as.ppp(longleaf)
a1 <- as.data.frame(longleaf)
#a1
#? as.ppp
#? as.psp
#? as.im
a

## Marked planar point pattern: 584 points
## marks are numeric, of storage type  'double'
## window: rectangle = [0, 200] x [0, 200] metres

plot(a)

aa <- as.ppp(chorley)
aa

## Marked planar point pattern: 1036 points
## Multitype, with levels = larynx, lung 
## window: polygonal boundary
## enclosing rectangle: [343.45, 366.45] x [410.41, 431.79] km

plot(aa)

#aa$window$bdry

knitr::include_graphics('scale1to10.jpg')

library(spatstat)
murchison

## List of spatial objects
## 
## gold:
## Planar point pattern: 255 points
## window: rectangle = [352782.9, 682589.6] x [6699742, 7101484] metres
## 
## faults:
## planar line segment pattern: 3252 line segments
## window: rectangle = [352782.9, 682589.6] x [6699742, 7101484] metres
## 
## greenstone:
## window: polygonal boundary
## enclosing rectangle: [352782.9, 681699.6] x [6706467, 7100804] metres

mpp <- as.ppp(murchison$gold)
mpp

## Planar point pattern: 255 points
## window: rectangle = [352782.9, 682589.6] x [6699742, 7101484] metres

plot(mpp)

lpp <- as.psp(murchison$faults)
lpp

## planar line segment pattern: 3252 line segments
## window: rectangle = [352782.9, 682589.6] x [6699742, 7101484] metres

plot(lpp)

ppp <- as.polygonal(murchison$greenstone)
ppp

## window: polygonal boundary
## enclosing rectangle: [352782.9, 681699.6] x [6706467, 7100804] metres

plot(ppp)

Moran’s I:

For continuous attribute \(Z\) with \(E[Z(s)] = \mu\) as well as constant variance:

\[I = \frac{n}{(n-1)s^2 w_{..}} \sum_{i=1}^n \sum_{j=1}^n w_{ij} (Z(s_i) - \bar{Z})(Z(s_j)-\bar{Z}), \textrm{ where } w_{..} = \sum_i \sum_j w_{ij}\]

• If \(I > \frac{-1}{n-1}\): a location tends to be connected to locations with similar attribute values.

• If \(I < \frac{-1}{n-1}\): attribute values at locations connected to a particular location tend to be dissimilar

• Local Moran (Anselin 1995): \[I_i = n (Z_i - \bar{Z}) \frac{\sum_j w_{ij} (Z_j - \bar{Z})}{\sum_i (Z_i - \bar{Z})^2}\]

library(raster)

## Loading required package: sp

## 
## Attaching package: 'raster'

## The following objects are masked from 'package:spatstat':
## 
##     area, rotate, shift

## The following object is masked from 'package:nlme':
## 
##     getData

#? Moran
r <- raster(nrows=10, ncols=10)
r[] <- 1:ncell(r)
image(r)

Moran(r)

## [1] 0.8353565

m <- MoranLocal(r)
plot(m)

com <- -1/(99)
com

## [1] -0.01010101

# Rook's case
f <- matrix(c(0,1,0,1,0,1,0,1,0), nrow=3)
Moran(r, f)

## [1] 0.8929367

Geary’s C:

For continuous attribute \(Z\) with \(E[Z(s)] = \mu\) as well as constant variance we have a type of autocorrelation measure:

\[C = \frac{1}{2S^2 w_{..}} \sum_{i=1}^n \sum_{j=1}^n w_{ij} \left(Z(s_i)-Z(s_j) \right)^2\]

If \(C>1\) the locations are connected to locations with dissimilar values and vice versa for \(C<1\).

The constant mean and variance is important: if this is not true the similarity/dissimilarity is more likely due to the heterogeneous mean and variance.

#? Geary
r <- raster(nrows=10, ncols=10)
r[] <- 1:ncell(r)
image(r)

Geary(r)

## [1] 0.04739631

plot(GearyLocal(r))

# Rook's case
f <- matrix(c(0,1,0,1,0,1,0,1,0), nrow=3)
Geary(r, f)

## [1] 0.03095362

The Variogram/Semi-variogram (Matheron, 1963):

\(\gamma(h)\) is half the average squared difference between points \(s_1\) and \(s_2\) separated by distance \(h\).

\[2\gamma(h)=2\gamma(s_1,s_2) = E[(Z(s_1)-Z(s_2))^2]\]

Empirically calculated as:

\[\hat{\gamma}(h) = \frac{1}{2|N(h \pm \delta)|} \sum_{i,j \in N(h \pm \delta) }|Z(s_1)-Z(s_2)|^2\] where \(N(h \pm \delta) = \{(s_i,s_j): |s_i,s_j| \in h \pm \delta, i,j = 1,2,..., N\}\).

YouTube Video: https://www.youtube.com/watch?v=SJLDlasDLEU&t=18s

• At least 30 pairs of locations should be available at each \(h\), thus group into \(h \pm \epsilon\).

library(sp)
library(gstat)

## Registered S3 method overwritten by 'xts':
##   method     from
##   as.zoo.xts zoo

## 
## Attaching package: 'gstat'

## The following object is masked from 'package:spatstat':
## 
##     idw

#? variogram
#? meuse
data(meuse)
#meuse$x
#meuse$y
plot(x = meuse$x, y = meuse$y)

# no trend:
coordinates(meuse) = ~x+y
vv1 <- variogram(log(zinc)~1, meuse)
plot(vv1)

# residual variogram w.r.t. a linear trend:
vv2 <- variogram(log(zinc)~x+y, meuse)
plot(vv2)

# directional variogram:
vv3 <- variogram(log(zinc)~x+y, meuse, alpha=c(0,45,90,135))
plot(vv3)

vv4 <- variogram(log(zinc)~1, meuse, width=90, cutoff=1300)
plot(vv4)

The semivariogram captures the second-moment structure of the data.

\[\gamma(s,h) = \frac{1}{2} [var(Z(s) - Z(s+h))] = \frac{1}{2} [var(Z(s)) + var(Z(s+h)) - 2cov(Z(s),Z(s+h))].\]

If second-order stationarity holds: \(\gamma(h) = \frac{1}{2} (2\sigma^2 -2 C(h)) = C(0) - C(h)\).

If positive spatial correlation is present:

• \(\gamma(0) = C(0) - C(0) = 0\). However due to measurement error and microscale variation there may be a nugget effect - this is often observed in sample variograms.
• The variogram attains \(\gamma(h^*) = \sigma^2\) when the lag exceeds the distance \(h^*\) for which the data points are correlated.
• The range is \(h^*\) (where \(C(h^*) = 0\) or at \(95\% \sigma^2\))

Empirical semivariogram:

\[\hat{\gamma}(h) = \frac{1}{2 |N(h)|} \sum_{N(h)} \left(Z(s_i) - Z(s_j) \right)^2 \]

Empirical covariance estimation:

\[\hat{C}(h) = \frac{1}{|N(h)|} \sum_{N(h)} (Z(s_i)-\bar{Z})(Z(s_j)-\bar{Z]})\]

knitr::include_graphics('variogram.jpg')

Autocorrelation function for a time series

For a time series \(Z(t_1),...,Z(t_n)\) with \(E[Z(t_i)] = 0\) and \(var(Z(t_i)) = \sigma^2, i = 1,...,n\) the covariance function is a function of the distance in time \(t_j - t_i\) \[cov(Z(t_i),Z(t_j)) = E[Z(t_i)Z(t_j)] = C(t_j-t_i).\]

The autocorrelation function is \[R(t_j-t_i) = \frac{cov(Z(t_i),Z(t_j))}{\sqrt{var(Z(t_i))var(Z(t_j))}} = \frac{C(t_j- t_i)}{C(0)}.\]

But time has direction, space does not…

Autocorrelation function in space

Covariance function: \(C(s,h) = cov(Z(s),Z(s+h)) = E[(Z(s) - \mu(s))(Z(s+h) - \mu(s+h))]\)

Autocorrelation function (correlogram): \[R(s,h) = \frac{C(s,h)}{\sqrt{var(Z(s)) var(Z(s+h))}}\]

Edge effects are more important than in time series since potentially many points located at the domain boundary.

Spatial dependencies may develop differently in various directions - anisotropy.

Second-order stationarity: \(Z(s)\) has a constant mean \(\mu\) and a covariance function which is a function of the spatial separation between points only i.e. \(C(h)\). (\(C(0) = cov(Z(s), Z(s)+0) = var(Z(s)) = \sigma^2\) so not a function of spatial location.)

So under second-order stationarity: \(cov(Z(s),Z(s+h)) = C(h)\) abd \(R(h) = \frac{C(h)}{\sigma^2}.\)

plot(murchison$gold)

den <- density(murchison$gold)
plot(den)

den2 <- density(murchison$faults)
plot(den2)

#? bw.diggle
#? bw.scott
#? bw.stoyan
#? bw.frac

den <- density(murchison$gold, bw.diggle)
plot(den)

den <- density(murchison$gold, bw.scott)
plot(den)

den <- density(murchison$gold, bw.stoyan)
plot(den)

den <- density(murchison$gold, bw.frac)
plot(den)

Semivariogram and Covariogram

Semivariogram: \(\gamma(h)\) vs. \(h\)

Variogram: \(2\gamma(h)\) (but naming is not consistant)

Covariance function: \(C(h)\) vs. \(h\)

Under stationarity \(\gamma(h) = \gamma(s_i-s_j) = C(0) - C(s_i-s_j)\). Recall that \(C(0) = var(Z(s))\).

ESTIMATION (Matheron 1962, 1963): \[\gamma(s_i-s_j) = \frac{1}{2|N(s_i-s_j)|} \sum_{N(s_i-s_j)} ({Z(s_i) - Z(s_j)})^2\]

\[C(s_i-s_j) = \frac{1}{2|N(s_i-s_j)|} \sum_{N(s_i-s_j)} (Z(s_i) - \bar{Z})(Z(s_j)-\bar{Z}) \textrm{ where } \bar{Z} = \frac{1}{n} \sum_{i=1}^n Z(s_i). \]